#pragma once

#include "Tensor.h"
#include "Polynomial.h"

/**
 * @file Tensor.h
 * @author xingyifan
 * @date 2022-07-26 18:26
 *
 * @description: 张量对应网格格点
 */

template <int Dim>
using Grid = Tensor<Dim>;

/**
 * @file Calculus.h
 * @author xingyifan
 * @date 2022-07-26 00:12
 *
 * @description: 被离散化的方形区域
 */

template <int Dim>
class Domain
{
public:
    Domain(Real *l, Real width) : width_(width)
    {
        for (int i = 0; i < Dim; i++)
        {
            lower_[i] = l[i];
        }
    }
    Domain(std::initializer_list<Real> l, Real width) : width_(width)
    {
        auto it = l.begin();
        for (int i = 0; i < Dim; i++)
        {
            lower_[i] = *it++;
        }
    }
    Domain(const Domain &other) : width_(other.width_)
    {
        for (int i = 0; i < Dim; i++)
        {
            lower_[i] = other.lower_[i];
        }
    }

    Real width() const { return width_; }
    int dim() const { return Dim; }
    const Real *lower() const { return lower_; }

    /**
     * @file Calculus.h
     * @author xingyifan
     * @date 2022-07-30 14:58
     *
     * @description: 计算区域体积
     */

    Real volume() const { return pow(width_, Dim); }

    /**
     * @file Calculus.h
     * @author xingyifan
     * @date 2022-07-26 00:13
     *
     * @description: t 当前指标对应的格点通过 Domain 转换为区域中对应的格点
     */

    Vector gridPoint(const Grid<Dim> &t) const
    {
        Real h = width_ / (t.base() - 1);

        // 获得角点
        Vector point(Dim, 0);
        for (int i = 0; i < Dim; ++i)
        {
            point[i] = lower_[i] + t.index(i) * h;
        }
        return point;
    }

    friend std::ostream &operator<<(std::ostream &out, const Domain &d)
    {
        for (int i = 0; i < Dim; i++)
        {
            out << "[" << d.lower_[i] << "," << d.lower_[i] + d.width_ << "]" << std::endl;
        }
        out << std::endl;
        return out;
    }

private:
    const Real width_; // 宽度
    Real lower_[Dim];  // 下界
};

/**
 * @file Calculus.h
 * @author xingyifan
 * @date 2022-07-26 00:12
 *
 * @description: 导数矩阵：FuncX 接收 Dim 元向量，输出 Dim 元向量
 */

class Jacobian
{
public:
    // 多元向量微分，返回 Jacobi 矩阵
    Matrix operator()(FuncX f, const Vector &x)
    {
        int n = x.width();
        Real h = 1e-8;

        Matrix d(n, 0);
        for (int i = 0; i < n; i++)
        {
            Vector y1 = x, y2 = x;
            y1[i] += h, y2[i] -= h;

            // 中心差分向量
            Vector dp = (f(y1) - f(y2)) / (2 * h);
            for (int j = 0; j < n; j++)
            {
                d[j * n + i] = dp[j];
            }
        }
        return d;
    }
};

/**
 * @file Calculus.h
 * @author xingyifan
 * @date 2022-07-26 00:13
 *
 * @description: 多元数值积分 D 阶，Dim 维
 */

template <int D, int Dim>
class Integrator
{
public:
    /**
     * @file Calculus.h
     * @author xingyifan
     * @date 2022-07-28 00:18
     *
     * @description: 在构造函数中初始化权重
     */

    Integrator()
    {
        // 计算权的分布
        Tensor<1> baseWeight(D + 1, 0);

        // 使用插值多项式积分计算权
        std::vector<Real> x(D + 1);
        for (int i = 0; i < D + 1; i++)
        {
            x[i] = i;
        }
        for (int i = 0; i < D + 1; i++)
        {
            std::vector<Real> f(D + 1, 0);
            f[i] = 1;
            baseWeight[i] = Polynomial::interpolate(f, x).intergral(0, D);
        }

        // 根据维数取权张量的积
        weight_ = TensorOp::power<Dim>(baseWeight);
    }

public:
    /**
     * @file Calculus.h
     * @author xingyifan
     * @date 2022-07-28 00:11
     *
     * @description: 向量函数积分：FuncX 接收 Dim 元向量，输出 n 元向量
     */

    Vector NewtonCotes(Domain<Dim> &domain, FuncX f)
    {
        // 分割宽度
        Real h = domain.width() / D;

        // 离散 Dim 维区域为网格
        Grid<Dim> grid(D + 1, 0);

        // 获得输出向量的长度
        int n = f(Vector(Dim, 0)).size();
        Vector I(n, 0);

        // 逐分量积分
        for (int i = 0; i < n; i++)
        {
            // 填充区域
            do
            {
                Vector y = f(domain.gridPoint(grid));
                // 获得对应区域中的点计算函数值
                grid.value() = y[i];
            } while (grid++);

            // 获得单元值并加权平均
            Tensor<Dim> value = grid.cell();
            I[i] = TensorOp::dot(weight_, grid) * pow(h, Dim);
        }
        return I;
    }

public:
    Tensor<Dim> weight_; // 积分权
};